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Abstract 


A numerical study of enhanced oil recovery from porous formations is presented in this work 
Resistance to oil movement arises from viscous as well as surface tension forces The goal 
of enhanced oil recovery techniques is to reduce these forces and improve the mobility of oil 
Thermal methods are frequently used to lower the oil viscosity, while surfactants are used to 
lower the surface tension The present study concerns modeling of the oil recovery by hot and 
cold water injection and the effect of surfactants on the enhancement of the oil recovery 

Equations governing nomsothermal two-phase flow in a porous medium are obtained 
from the mass balance principle for the two phases, Darcy's law as applicable to an unsatur- 
ated porous region and the energy equation in terms of a local volume-averaged temperature 
These are supplemented by constitutive relations that connect the capillary pressure, relative 
permeability and the water saturation The governing equations are reduced to a coupled 
system of nonlinear partial differential equations in terms of oil and water pressures and tem- 
perature In the present study, pressure equations are solved using an implicit central difference 
scheme and saturation is determined from the difference between oil and water pressures The 
energy equation has been solved using an operator splitting (OS) algorithm 

A domain decomposition technique has also been developed in the present work for 
solving the governing equations Domain decomposition is a useful method for solving problems 
over very large domains or regions having a complex geometry An attractive feature of this 
technique is the possibility of parallelizing the computer code suitable for computers with a 
parallel architecture In the present work, a completely parallel code has been developed for 
the oil recovery problem using the domain decomposition technique Results do show a great 
potential for this technique for simulating the field scale problems 

The present study also concerns modeling of flow in a thick porous layer considered as 
a two-spot model to predict the actual flow behaviour that takes place in a five-spot model 
Results reveal the flow pattern of the field scale problems 
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Chapter 1 


Introduction 


It has been recognised that the mankind is likely to face a critical shortage of petroleum 
all over the world in near future Under such a situation, exploring the possibility of 
enhancing the oil recovery from the existing reservoirs invariably assumes tremendous 
importance It is also well known that a large proportion of the existing oil reservoirs are 
still underrecovered The need for increasing the recovery of oil from such reservoirs calls 
for developing the enhanced oil recovery (EOR) techniques Particularly m our country, 
the availability of oil will play a vital role in the nation’s competetiveness and industrial 
health over the next two decades Several natural resources of oil have been identified in 
the country, but all of these are being operated at the primary stage, with no enhancement 
technology Clearly there is a need for establishing the technology of EOR at the earliest 
And the objective of the present work is to consider the EOR techniques from numerical 
point of view and develop faster and more reliable numerical codes for reservoir simulation 

The conventional methods of oil recovery, earlier to the advent of EOR techniques, 
can be categorised as the following 

(1) Primary recovery In this case the m-situ pressure m the reservoir is often high enough 
to force the resident fluid out of the production wells without any pumping effort This 
accounts for 15%— 20% of the total oil reserve 

(2) Secondary recovery Here a fluid such as water is injected into some wells m the reser- 
voir while petroleum is pushed through the production wells This serves the dual purpose 
of maintaining high reservoir pressure and of flooding the porous medium to physically 
displace some of the oil and push it towards the production wells But unfortunately, 
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waterflooding is still not extremely effective and significant amounts, upto 50% of hydro- 
carbon often remain in the reservoir Due to strong surface tension effects, a large amount 
of oil is trapped m small pores with narrow throats and is not washed out with routine 
waterflooding techniques 

In order to recover more hydrocarbon, several enhanced oil recovery techniques in- 
volving complex chemical and thermal effects have been developed A brief overview of 
various EOR techniques is furnished here 

(1) Thermal recovery methods (Boberg, 1988, Aziz, 1986) aim at reducing the oil 
viscosity by raising its temperature, thus improving the displacement efficiency of oil It 
involves injection of either hot water or steam into the reservoir via injection wells As 
the heated fluid comes into contact with the cooler oil, the temperature of oil is raised 
Consequently the oil viscosity is reduced and it is pushed rather easily through the porous 
formation Pushing the high temperature fluid down into the reservoir is a difficult engin- 
eering problem The hot injected fluid may lose heat while moving through the injection 
well and later to the host rock bounding the oil rich region Some energy is also consumed 
m the chemical reactions of hydrocarbons Apart from heat loss, one major problem is 
the gravity override This occurs when the hot fluid rises on to the top of the oil saturated 
sand layer and bypasses it If the flow rate is high, the interface between the resident 
petroleum and the invading displacing fluid can be unstable leading to the formation of 
long fingers These fingers grow m length towards the production wells, thus bypassing 
much of the oil Such factors cause a drastic reduction m oil production 

(2) Surfactant flooding improves the oil displacement efficiency by reducing the mterfa- 
cial tension between oil and water considerably Surfactants are surface active chemicals, 
a solution of which m water, injected at an appropriate pressure mobilizes the trapped 
oil to form a flowing oil bank by overcoming the capillary forces that are predominant m 
unsaturated porous regions 

(3) Polymer flooding (Boberg, 1988, Shah et al , 1977) is used to increase the volume of 
the reservoir contacted by the displacing fluid Addition of a low concentration of polymer 
makes the floodwater more viscous, decreases its tendency to finger through and bypass 
oil, and leads to increased sweep of the reservoir. Polymer flooding, however, does not 
recover residual oil trapped in microscopic pores by capillary forces Therefore, it is very 
often viewed as an improved secondary recovery process rather than an EOR technique 

(4) In-situ combustion technique is generally employed for reservoirs that are too deep 
to be mined and it calls for m-situ transformation of the hydrocarbons into states which can 
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be pumped to the surface via production wells. In this technique, the solid hydrocarbon 
is ignited m the ground and oxygen or air is pumped into the injection wells to maintain 
combustion The hot hydrocarbon gases are pumped to the surface via the production 
wells for use as low-grade hydrocarbon 

Reservoir simulation is one of the most challenging problems m engineering (Ewing, 
1983, Jim Douglas Jr , 1983). Owing to variety of contributing factors like the inherent 
complexity of the domain, complex nature of the constitutive relationships and instability of 
the oil-water interface the governing partial differential equations are highly non-linear and 
coupled Also analytical solutions are difficult to obtain. In the present work, we propose 
to perform the numerical simulation of the EOR technique using hot water injection. This 
aspect of the investigation is an extension of the earlier work (Pillar and Muralidhar, 
1993) We also propose to investigate the effect of surfactants on oil recovery for a two- 
dimensional rectangular porous formation This part of the investigation is in continuation 
with the earlier work of Chatterjee and Muralidhar (1995) Also a numerical study has 
been done for a two-dimensional two-spot model, more precisely known as quarter five- 
spot model (Sorbie et al , 1995, Zhang et al , 1997) For the two-spot problem, we have 
studied the isothermal injection and also the effect of surfactants on oil recovery. For the 
numerical investigations mentioned above, a pressure based formulation (Ewing, 1983) has 
been used In the present study, pressure equations are solved using an implicit central 
difference scheme and saturation is determined from the difference between oil and water 
pressures The energy equation has been solved using an operator-splitting algorithm 
(Muralidhar et al , 1993) The presence of mass transfer phenomenon such as adsorption 
and desoiption of surfactants from water to the solid matrix and vice versa is neglected in 
the present work 

Computational time is an important issue while solving reservoir problems Because 
of high degree of non-linearity of the governing equations and large size of the reservoirs, 
reservoir simulation is computationally intensive. However, the current approach towards 
handling such complex problems is to improve the numerical algorithm on one hand and 
the computer performance is on the other With the advent of parallel computers there 
exists a great potential in increasing computer performance m terms of speed as well as 
memory In this context, domain decomposition techniques (Glowmski et al , 1983 , Le 
Tallec, 1994) are very suitable for reservoir simulation Splitting the physical domain 
into sub-domains is the essential feature of the domain decomposition. The governing 



4 


equations are then solved independently m each sub-domam The individual solutions are 
subsequently assembled carefully to get a converged solution of the original problem on 
the complete physical domain Domain decomposition can rapidly make the simulation 
parallelizable since calculations m each sub-domam can be carried out independently on 
a processor They can also handle large problems on sequential machines more efficiently 

A great amount of investigations on domain decomposition have appeared in the 
literature m the past few years Glowmski et al. (1983), have developed a number of 
efficient parallelizable domain decomposition algorithms m their pioneering work and these 
are m wide use Perng and Street, (1990) have applied domain decomposition technique 
to solve Navier-Stokes equations in complex domains which are otherwise difficult to solve 
without coordinate transformation Parallelizable domain decomposition codes applied 
to transient problems are however scant in literature Sequential domain decomposition 
codes applied to steady problems are available in plenty (Perng et al , 1990) 

In the present work, we have developed a parallel code to solve the transient oil 
recovery problem for isothermal injection and the effect of surfactants (Muralidhar et al , 
1996) has also been studied Uzawa’s algorithm (Glowinski, 1983, Yagawa, 1991) has been 
used for interface treatment m conjunction with an implicit formulation The algorithm 
has the desired parallel and temporal properties Results are in excellent agreement with 
those due to the full domain simulation 

The present work is organised along the following sequence. 

• Mathematical modelling and formulation of the oil recovery problem. 

• Overview of the numerical techniques developed and implemented for the simulation 

• Results and discussions of the simulation of the EOR problem 

• Conclusions and scope for the future work. 



Chapter 2 


Mathematical Model and Formulation 
of the Problem 


In this chapter the mathematical formulation of the governing equations for two-phase flow 
through a porous media has been discussed The physical laws that govern multi-phase 
flow through the porous medium are 

1 Conservation of mass of the individual phases, 

2 Appropriate form of Darcy’s law that includes capillary effects, and 

3 Conservation of thermal energy 


2.1 As sumptions 


The following assumptions have been made in the present study 

1 The fluid phases, namely oil and water flow simultaneously 

2 Since the phase velocities are small, local thermal equilibrium between the fluid and 
solid phases is assumed 

3 Gravity term is neglected m the momentum equation It is reasonable because the 
time frames axe short and as a result the gravity override effect will not show up 

4 The fluids, the solid matrix and the porous region are assumed to be incompressible 
and porosity is a constant 

5 The fluids are immiscible, i e , there is no mass transfer between them 
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6 Relative permeabilities are taken as unique functions of saturation and viscosities as 
unique functions of temperature 

7 Specific heats are assumed to be constant 

8 There is no net mass source m the flow field 

9 Viscous dissipation and mechanical work terms in the energy equation are neglected 

10 Fluids, l e oil and water are assumed to exhibit Newtonian behaviour 


2.2 Governing Equations 

Equations governing transient two-phase flow through a homogeneous, isotropic porous 
medium are given as follows 

Mass Balance 

Accumulation rate + Efflux rate = Generation rate 
i e 

d 

eS t p ,) + V (, p t u t ) = 0 (l = oil, water) (2 1) 

(assuming no generation) 

Darcy’s Law 

The momentum conservation equations for low Reynolds number, i e Re<l, assume 
the form of Darcy’s laws which are given as follows 

«. = Vp, (2 2) 

llx 

(neglecting the effect of gravity override) 

Energy Equation 

The rate at which thermal energy diffuses into a control volume is equal to the sum 
of the rate of accumulation and its efflux rate 
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^- + UT X + VT y = — V 2 T 
at ax 

where 

jj UopoCo ^'wPvjC'w 
&T 

' 1^. ^oPo^o "h V w p w C w 

(Tt 

(JX = ^{(1 *h + (1 c) (p^)# 


(2 3) 


Constitutive relationships 


^ _ 2_ / flgA 

p» \ dp* ) J’ 




(2 4) 
(2 5) 


(from the definitions of the coefficient of compressibility and the coefficient of expansivity 
respectively) 

Pcow (Syj) = Po Pw (2 6) 


(p Cow is capillary pressure and a unique fimction of S w ) 


S 0 + S w = 1 (2 7) 

(Here we assume that the pore volume not occupied by oil, is occupied by water and no 
gas phase is present ) 


k rt — k rt {S-uj'j 

(l e relative permeability of a phase is a unique function of saturation ) 


(2 8 ) 


Equations 2 1 and 2 2 can be combined to generate two equations for oil and water 
respectively, using the constitutive relationships These equations are as follows 


— SB — 4 - S £ ^ Pw 
boP° dt + B o l; 0 dt 


■S 0 — + S £ ^ Pw + 

wHw i i 


dS w 

dp 0 

dp w ‘ 

= Iv 

K kroPo 

dPc ow 

. dt 

dt m 

Po 



dS w 

dp 0 

& 
_1 

_J_ V 

K k rw p w ^ 

dPcow 

. dt 

dt . 

Pw 

■ 

vo 

$ 

=1 

1 


^ Pw 


(2 9) 

(2 10 ) 
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It is extremely important to justify the use of the pressure formulation at this stage 
Oil recovery problem can also be formulated in terms of water saturation This means that 
the three dependent variables which will appear m the coupled system of non-linear partial 
differential equations will be in this case oil pressure, water saturation and temperature 
instead of oil and water pressures and temperature However this approach can lead 
to difficulties m numerical computation since the oil-water interface can be quite sharp 
In contrast to this, the pressure field must be continuous everywhere This results in 
appreciable improvement m the accuracy of the numerical simulation 


2.3 Physical Description 


The problem of interest considers flow through a thm porous layer and a thick porous 
layer 


The physical domain for the flow through a thm porous layer is shown in Figure 2 1 
The flow is two dimensional and the rectangular domain has a nominal size of 10 m x 
12 m. Essentially we are studying the water piston effect in a rectangular domain packed 
with sand and guarded with impermeable but thermally conductive rocks above and below 
the porous region (called overhang and underhang respectively) 

The physical domain for the flow through a thick porous layer is shown in Figure 2.2 
The geometric model is the much talked about five-spot model (Figure 2 2 (a)) The flow 
is essentially two-dimensional and the square domain (the plan view of the thick porous 
layer) has a nominal size of 6 m x 6 m We have, however, studied the flow for the 
quarter five-spot model (hereafter mentioned as two-spot model; Figure 2 2 (b)) We aim 
at simulating the flow for a set of three simplified versions (Figure 2 2 (c), (d) and (e)) 
of the original two-spot model Barring the regions for injection and production wells, the 
sand packed domain is guarded with impermeable rocks 

Initially the porous medium is impregnated with oil to a given saturation It has 
some initial temperature, called the formation temperature and is at rest One of the 
assumptions in this study is that the pore volume not containing oil has water. At equi- 
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librium, the oil and water pressures in the porous layer can be uniquely determined from 
their saturation Oil recovery is initiated by pushing water at a given higher pressure and 
a given temperature We observe the movement of saturation and temperature profiles m 
the porous region 



p„ = 1 31 MPa 


If = 50°C 

(Oil + Water) 
Extraction 


Figure 2 1 Description of the Oil Recovery Problem from Thm Porous Formation 
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Figure 2 2. Description of the Oil Recovery Problem from Thick Porous Formation 
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We have addressed to two variants of the mam physical phenomena m this study as 
furnished below 

(1) Isothermal Problem In this case the high pressure water is injected at a temperature 
same as that of the formation temperature The effect of surfactants is also considered for 
this case We have accomplished this study for both the thm porous formation and the 
thick porous formation Essentially we need not solve the energy equation herein 

(2) Non-isothermal Problem The problem of interest involves transportation of thermal 
energy into the flow domain through hot water that is injected at a higher pressure Nat- 
urally the temperature distribution is important in this case and the energy equation is 
solved along with the pressure equations However, we have undertaken this study only 
for the thin porous formation 


2.4 Initial and Boundary Conditions 


The pressure equations (equations 2 9 and 2.10) are solved for both the rectangular and the 
square flow domains subject to a quiescent initial condition, impermeable confining walls 
and specified oil and water pressures at the inflow and outflow sections At the injection 
zone of the domain, the water is at a pressure of 1 7926 MPa and at a temperature of 
100° C Saturation of water is 0 86 here It is not unity because of the presence of some re- 
sidual oil which cannot be flushed out It is by far the maximum water saturation possible 
anywhere in the domain. At the production zone of the domain, water pressure, temper- 
ature and saturation are maintained at 1.31 MPa, 50° C and 0 2 respectively Obviously 
the injection water temperature is 100° C only for the non-isothermal injection, otherwise 
the temperature is same as that of the production zone 

Initially, the porous region is maintained at a uniform water pressure and temperature 
of 1 31 MPa and 50°C respectively and the oil saturation is of 0.8 Oil pressure is obtained 
by adding water pressure to the capillary pressure, the latter being a unique function of 
the water saturation. 

The impermebility condition at the bounding surfaces calls for the Darcy velocity 
components becoming zero, which give rise to the boundary conditions as, 
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dp, dp, n 

8x 8y 

Boundary condition related to heat loss to overhang and under hang correpondmg to the 
thin porous layer flow domain is given by 

csrp 

-^j + Bi(T-T f ) = 0 at y = 0, Y 

The model presented in this work, predicts the movement of the water front from 
the injection zone into the oil-rich porous region This increase m water saturation is a 
measure of the amount of oil displaced Oil displacement is measured m terms of the 
percentage pore volume (ppv) It is defined as, 

, total volume of oil produced 

ppv = 100 x — — 

total pore volume 

The use of Dirichlet boundary condition at the outflow section is justified only if the 
location of this section is far away from the injection region In the present study, oil 
recovery is computed over half a distance of the domain size For a five hour simulation, 
considered here, numerical experiments show that oil recovery is insensitive to the location 
of the outflow section 

The data for the parametric values and constitutive relationships used m the present 
study have been mostly adapted from Boberg (1988) and are given m Appendix A 



Chapter 3 


Numerical Scheme 


In this chapter we shall discuss, in detail, the numerical techniques used for solving the 
governing differential equations enumerated earlier Also the implementation of the domain 
decomposition technique will be taken up in our discussion 


3.1 Discretization of Pressure Equations 


The oil and water pressure equations 2 9 and 2 10 are discretized using a control volume 
finite difference scheme for both the flow domains (Figure 2 1 (a) and Figure 2 2 (b)). The 
scheme is fully implicit m time This implies that all spatial derivatives are evaluated at the 
new time step The time derivatives of pressure are discretized using the forward difference 
technique and those of temperature using the backward difference technique However the 
backward differencing in time of the temperature term is justified because the effect of 
temperature on the pressure distribution is brought about through the fluid properties 
The fluid properties do not change as rapidly with time as the primary variables like phase 
pressures and saturation The final form of the discretized pressure equations are 

[Aopo t+h , + B oPotJ + C 0 p Wx] + D 0 p 0t _ l} + E 0 p 0uj+l + F 0 p 0tJ _ J n+1 = G 0 (3.1) 


\A w p Wt + T B w p Wtj T C w p 0t j T D w p lUi _ l j T E w p WlJ _^ 1 T F w Pw t ^i — G w (3-2) 
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hJ 
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where 


K kropo 
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K krwPw 

Cfinj 


Here (9 e ,9 w ,9 n ,9 s ) and (7 e ,7u„7n,7s) are values corresponding to the points lying 
towards the east, west, north and south of the node i,j (as shown m the Figure 2 1 (b)) 
and are evaluated as the harmonic averages of the values at nodes that he on either side 

The oil pressure equation 3 1 and the water pressure equation 3 2 are written for 
every node on the grid mesh to obtain a set of 2 x n x m simultaneous eqations These 
equations can be written m a matrix form The matrix thus obtained is banded The 
matrix structure is shown m Figure 3 1 It has been solved by using a specialized sparse 
matrix solver based on Gaussian elimination (Duff, 1980) and also by preconditioned 
conjugate gradient method (PCG; Muralidhar and Sundararajan, 1995) 
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Figure 3 1 Matrix Structure in Pressure Calculation 


3.2 Solution of the Energy Equation 


The energy equation is solved by using the operator splitting algorithm, The energy 
equation, 



dt d X dy Ur/ 

(3 3) 

is split into the following steps 


Predictor: 

dT rr dT ,,8T „ 

aF + £/ a? + y 6fe-° 

(3 4) 

Corrector: 


(3 5) 


3.2.1 Solution of the Predictor Step Using Streamlines 

We choose a curvilinear coordinate system (£, 77 ) locally at each node such that V£ is 
aligned with the net fluid composite velocity at that node, i e £ is streamline passing 
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through the node (Figure 3 2 (a)) Then 

UT X + VT y = U% (3 6) 

where, U' = VU 2 + V 2 

Hence the predictor equation becomes 

T t + U% = 0 (3 7) 

The solution of this equation is T(U't — £) = Constant 

The solution of the Predictor step requires that the temperature at Q at the n th timestep 
be transferred to the point P at (n + l) th time step (Figure 3 2 (b)), l e 

Tp\U\t + At) - £) = T$(U't -((- Af)) 

where, A£ = U ,n A t 

The temperature Tq is obtained through interpolation 
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Figure 3 2 Solution of Predictor of the Energy Equation. 



3 2 Solution of the Energy Equation 


18 


3.2.2 Solution of the Corrector Step Using ADI 


The corrector step consists of the two-dimensional unsteady state heat conduction equa- 
tion It is solved using the Alternating Direction Implicit (ADI) technique In ADI, 
we first compute row by row intermediate temperatures while sweeping along the j-axis 
Then by sweeping along the z-axis final temperatures for the time step are calculated in 
the column by column fashion The time step for each sweep is (A t) a i % = At/2 The 
discretization is given as 


Sweeping along j-axis: 


rp n +2 

J. r 



(A t)gdt 

Pe 




+ T t _ 


1,3 


2 T t 




(Ax) 7 


n+§ 
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(Ay) 2 


(3 8) 


or 


{AT t „ hj + B left Ti, j + CT t+hj } n+ > = {B rtght T hJ + DT hJ+1 + ET t ,. x } n (3 9) 
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{AT,,,-,. + B h „T,,, + CT„ J+1 }” +I 


= {Br„uT tJ + DT,+i,, + 


(3 11) 


where 


A = C = 


(At) adt 

Pe(Ay) 2 ’ 


Bleft — 1 + 


2(At)ad.i 

Pe(Ay) 2 ’ 


r-) — P — (A t)adt o _ 1 2(Af)(id, 

~ p e ( Ax) 2 ’ “ 1 Pe(Ax) 2 

The equations are written for every node m the row or column and the set of simul- 
taneous equations thus obtained yields a tridiagonal matrix, which can be inverted using 
TDMA as a matrix solver 


3.3 Algorithm 


The system of equations governing the distribution of oil and water pressures and tem- 
perature are non-linear and mutually coupled They are solved simultaneously and by 
iteration 

The algorithm, m brief, is as folows 

(1) The initial distributions of p oy p w , T and S w at t = 0 are prescribed 

(2) The coefficients (A 0 — G 0 and A w — G w ) of the prssure equations are computed using 
the values of p 0 , p w , T and S w at the current time step and iteration level 

(3) The system of pressure equations is solved to obtain the new values of p 0 and p w at 
the nodes formed by the grid 

(4) S w , p 0 and p w are updated via the constitutive relations using the new pressure values 

(5) Darcy velocities are calculated using the latest values of p 0 , p w , S w , p 0 and p w 

(6) Using S w , p 0 , p w and the Darcy velocities, the coefficients of the energy equation are 
computed and the equation is solved using the OS algorithm which yields new temperature 
distribution 

(7) p 0 and p w are updated once again using the new value of temperature The fluid 
viscosities are also updated at this stage 
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(8) Steps (2) — (7) are repeated until convergence of p 0 , p w and T is achieved 

(9) Fresh computation is initiated for the next time step starting from step (2) 

The above algorithm is used when the energy equation is involved, but for the iso- 
thermal injection, steps (5), (6) and (7) invariably become redundant 


3.4 Operator Splitting Algorithm 

3.4.1 Introduction 

The operator splitting algorithm (OS) is a special case of splitting methods or fractional 
step methods which reduces the solution of a complicated problem to a successive solution 
of simpler problems Various forms of splitting methos, such as, geometric splitting, phys- 
ical splitting and analytical splitting are available The OS algorithm used for the present 
investigation, belongs to the analytical splitting family This algorithm properly identifies 
the mixed mathematical character of the governing differntial equation and solves each ho- 
mogeneous component of the equation as accurately as possible Hence the OS algorithm 
solves the governing equation true to its character 


3.4.2 Model Example 

Let us consider a typical convection-diffusion problem (Muralidhar et al , 1993) arising m 
heat transfer problems of the form, 

^ + (u V)T = aW 2 T (3 12) 

The mathematical character of the equation is elliptic if t— > oo, parabolic if u is small and 
hyperbolic if u is large m magnitude In all other cases, the equation is said to be mixed 
m character These special cases are summarised below 

Steady state ( t -» oo) (u V)T = aV 2 T (elliptic) (3 13) 

Conduction limit (| u | — > 0) T t = aV 2 T (parabolic) (3 14) 

Convection limit (| u | -» oo) T t + u VT = 0 (hyperbolic) (3 15) 
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The OS algorithm handles this problem along the following line 

(1) At any time level t , equation (3 15) is solved over a time step At This is the predictor 
step 

(2) The equation (3 14) is then solved over the same time interval At and this is considered 
as the corrector step 

At high values of velocity, the correction required is quite small and and at low 
velocities the predictor step acting alone is inadequate However, the two steps combined 
will furnish the complete solution of the equation (3 12) irrespective of the magnitude of 
u 

3.4.3 Remarks 

The salient features of the OS algorithm are depicted below 

1. The scheme is completely free of upwmding and hence devoid of false diffusion errors 
that usually occur at high Peclet numbers 

2 The only error arising in the OS algorithm is due to time discretization and it is of 
the order At When either of the predictor or corrector equation is solved numeric- 
ally on a grid mesh of size Ax, then additional truncation errors (~ Ax n , n > 1) 
creep m Hence the overall error is, 

e ~ O (Ax n , At), n > 1 

where, n depends on the scheme used to solve the diffusive corrector step In our 
case for the corrector step, we have used a second order finite difference scheme and 
naturally the discretization error will be, 

e ~ O (Ax 2 , At) 

3. The main source of error in solving heat transfer problem lies in approximating the 
hyperbolic terms (u.VT) In the OS algorithm, these terms are treated analytically 
and solved as accurately as possible. In the oil recovery problem it has been solved 
on a local grid (Section 3 21) 
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4 On non-dimensionalisation, the heat transfer equation becomes 

dT ~ 1 

+ (u V)T = — v 1 2 r 

where, Pe is the Peclet number and is a measure of relative importance of convective 
transport with respect to diffusive transport It has been found from numerical 
experimentation that most of the existing algorithms fail or become inaccurate when 
Pe is high The OS algorithm has been found to be uniformly valid over any range 
of Pe As the Pe is raised, a discontinuity m the initial condition will propagate 
through the flow domain unchanged and is called a front Also fronts can be formed 
m non-lmear problems at low or moderate Peclet numbers The oil recovery problem 
considered here is one such problem 

5 The OS algorithm is conditionally stable since it involves explicit time marching and 
is limited by the Courant-Fnedrichs-Lewy (CFL) stability criterion 

3.5 Preconditioned Conjugate Gradient Method 

3.5.1 Introduction 

The Conjugate Gradient method can be used for inverting symmetric positive definite 
matrices It combines features of both iterative and direct solvers It converges faster 
than the Gauss-Siedel iteration for a given initial guess, but needs more arithmetic opera- 
tions per iteration This method minimises an error functional along mutually conjugate 
directions in each iteration Hence for a functional defined in N dimensions, it converges 
theoretically m N iterations For a matrix with a small spectral condition number, conver- 
gence is attained m just a few iterations, much less than N This method is particularly 
suitable for sparse matrices The CG algorithm stated for a matrix equation AX = b is 
as follows 


1 Assume X = X° (initial guess) 

2. Compute residue 

r k = b - Ax k 
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3 Set the direction vector 


pk = r k 


4 Compute steepest descent parameter, a 

( r k ,r k ) 

Q_ (P fc ,AP*) 

5 Update the values of x and r 

x k+1 = x k + ap k 
r k+i =r k_ aAp k 

6 Compute f3 


p 0V‘) 

7 Update the direction vector 

pk+l _ r k + 1 _|_ ppk 

8 Check for convergence in x 

9 Repeat steps 4 to 8 till convergence 


The above algorithm cannot be applied directly to a matrix which is not positive def- 
inite and symmetric Preconditioning is required for such matrices The rate of convergence 
depends upon the the spectral condition number of the matrix The Spectral condition 
number is defined as the ratio of maximum eigenvalue to the minimum eigenvalue of the 
matrix Closer it is to unity faster would be the convergence This requires the eigenvalues 
to cluster around a certain number. The aim of preconditioning is to improve the spectral 
condition number and move it towards unity. It involves the choice of a suitable positive 
definite matrix M so that AX = b may be replaced by M~ x AX = M _1 6, where M~ l A 
should have a large number of its eigenvalues closely grouped 

Three preconditioning strategies namely, (a) preconditioning I (symmetric matrices), (b) 
preconditioning II (asymmetric matrices) and (c) diagonal scaling have been described 
here 
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3.5.2 Preconditioning I (symmetric matrices) 

The equation AX = b is written as ( LL T )~ 1 AX = (LL r )~ 1 & where L is the lower triangu- 
lar matrix, which may be computed efficiently by Cholesky decomposition The Cholesky 
method is suitable only for symmetric matrices since an asymmetric matrix can not be 
factored as LL T If an approximate L matrix is used, Incomplete Cholesky Conjugate 
Gradient method (ICCG) is generated The complete algorithm is as follows 

V 

1 Assume X = X° (Initial guess) 

2 Compute residue 

r k = b - AX° 

3 Set the direction vector 

P*=(i T £)' , A T (lI T )' 1 r* 

4 Compute steepest descent parameter, a 

_ (r k ,(LL T )~ 1 r k ) 

( P k , (L T L)P k ) 

5 Update the values of x and r 

r k+ 1 = r k - aAP k 
X k+i = X k + aP k 

6 Compute /? 

(r*+‘, (££*•)-■,■*■») 

P (r k ,(LL T )~ 1 r k ) 

7 Update the direction vector 

P fc +! = [(L T L)- 1 A T (LL T )]- 1 r k+1 + (3P k 


8 Set k = k + 1, check for convergence m X. 

9 Repeat steps 4 to 8 till convergence 
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3.5.3 Preconditioning II (asymmetric matrices) 

For asymmetric non singular matrices the cholesky decomposition is replaced by the LU 
decomposition where we take unit diagonal elements in the upper triangular matrix U 
If L and U used have the same structure as A then the LU Decomposition Conjugate 
Gradient Algorithm (ILUCG) is generated Complete algorithm is as follows 

1 Assume X = X° (Initial guess) 

2 Compute residue 

r k = b - AX° 

3 Set the direction vector 

P k = (uTuy 1 A T (LL T )' 1 r k 

4 Compute steepest descent parameter^ 

__ (r k ,(LL T )~ 1 r k ) 
a ~ ( P k , ( U T U)P k ) 

5 Update the values of x and r 

; r k+1 =r k - aAP k 

X k+i =x k + aP k 

6 Compute /3 

( r k +\(LL T )- 1 r k+1 ) 

P ~ (r k , (LL T )~ l r k ) 

7 Update the direction vector 

pk+i = [(U T U )~ 1 A T (LL T )]- 1 r k + 1 + / 3P k 

8 Set k = k + 1, check for convergence m X. 

9 Repeat steps 4 to 8 till convergence 

3.5.4 Diagonal Scaling 

This strategy involves the scaling of matrix A so that the diagonal entries are strengthened 
with respect to the other elements within the matrix. 


D~ 1/2 AD~ 1/2 D lf2 X = D~ 1/2 b 
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AX = h 

A = D l/ 2 AD l/2 
X = D l/2 X 
b = D~ 1,2 b 

where D is the matrix containing only the diagonal entns of A 

We have applied PCG for solving the banded matrix arising out of the pressure 
equations as a substitute to the specialized sparse matrix inverter based on Gaussian 
elimination and the results are seen to be matching nicely 


3.6 Introduction to Domain Decomposition for 
Solving PDEs 

3.6.1 Opening Remarks 

Practical engineering problems inevitably require intensive computations because of the 
large size of physical domains and their complex geometries In many cases the com- 
putations cannot be implemented on existing computers owing to limitations of memory 
The current approach towards handling complex problems is to improve the numerical 
algorithm on one hand and the computer performance on the other With the advent of 
parallel computers there exists a great potential for increasing computer performance in 
terms of speed as well as memory However, parallel computers require parallel algorithms 
to make use of the distributed nature of the processor architecture Domain decomposition 
techniques are the simplest of parallel algorithms that can convert a sequential algorithm 
to one suitable for parallel computing Besides the ability to parallelize, other specific 
advantages of domain decomposition are 

1 It can generate smaller matrices that require less memory and are rapidly invertible 
when iterative techniques are used 

2 A complex assemblage of components of a system can be systematically analysed 
3. Complex phenomena can be localized by assigning separate subdomains to them 
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3.6.2 Domain Decomposition and Uzawa’s Algorithm 

The basic approach of domain decomposition is that a large domain is divided into many 
subdomains that are linked at the interfaces The governing equations are then solved 
independently in each subdomain with assumed interface conditions The efficiency and 
accuracy of domain decomposition lies in its interface treatment Among the various in- 
terface treatments available for domain decomposition, Uzawa ’s algorithm has been shown 
to be robust and parallelizable It is analytically well-behaved and shows monotone con- 
vergence properties 1 

Uzawa’s algorithm for an elliptic problem is presented below. Consider the linear 
governing differential equation 


A u — f m fi (3 16) 

and u = g on <9fl 

Here / and g are known functions, fl is the physical region, dVl is its boundary (Figure 3 3) 
and A is the Laplacian Operator Figure 3 3 also shows a decomposition of H into two 



Figure 3 3 Physical Region and its Division into Two Subdomains 

subdomains, 1 and 2, with C as their interface For a globally converged solution both u and 
its normal derivative are required to be continuous on C A start is made with a guess for 
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the derivatives, and (= A), where ni and n 2 are the normals at C extended into 
subdomains 1 and 2 respectively Using A^ as the boundary condition on C the governing 
differential equation is solved m regions 1 and 2 Let these solutions respectively be 
and U 2 ^ Let y\^ and j/ 2 ^ be the values of and u 2 W respectively on C Then A^ 1 ) 
is updated as 


A<« = A<‘> + 



(3 17) 


Here u> is an adjustable constant that lies between 0 and 1 and L is a suitable length scale 
It is taken to be the unity in the present work Using A^ 2 ^ as the next guess, the entire 
process is repeated After / iterations, Uzawa’s algorithm for updating the interface flux 


is 


A<™> = x m + 

L 

Interface convergence is judged by the criterion 

|a<w>-aM|<* 


(3 18) 


where 5 is the specified upper boimd The iterations of the algorithm given above converge 
when the computed function u^(x) matches the steady state solution of Equation (3 16) 

For transient problems, the approach adopted has been to apply the same algorithm 
between the successive time steps Each subproblem m the subdomain is solved by an 
implicit finite differnce scheme The choice of the implicit scheme guarantees the stability 
of the time-marching procedure and rate of convergence is dictated by the magnitude of the 
interface parameter p alone When the gradient of the dependent variable converges, the 
variable at the interface becomes continuous, and this is implied by equation (3.18) Hence 
the algorithm ensures continuity of u as well as its normal derivative at the interface. The 
only difficulty with Uzawa’s algorithm is in the choice of interface parameter p and must 
come from numerical experiments. 

Equation (3.18) shows that domain decomposition introduces interface iterations, 
even m linear problems Hence the possibility of an increase m CPU time must be con- 
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sidered, especially under transient conditions But it has been seen that actually the CPU 
time decreases The reason for obtaining this result is the following When matrix in- 
version is accomplished using iterative methods, for example the Gauss-Siedel method, 
the convergence rate varies as N ~ a , where N is the matrix size and a is some positive 
constant Hence convergence is significantly delayed for large values of N Application of 
domain decomposition cuts the matrix size by a factor, >2, depending on the number of 
subdomains used Hence each subdomain calculation is accelerated, providing an overall 
computing time advantage 


3.6.3 Domain Decomposition for Non-linear Problems 

In case of a non-lmear problem, the interface treatment must be modified Let us take 
up an example to illustrate this point Consider the linear heat conduction problem in the 
same domain as in Figure 


KV'T = P c p — (3.19) 

where K, p and c p are constants In Equation (3 19) the boundary conditions to be 
matched at the interfaces of the subdomains would be the temperature and the temperature 

0/71 

derivative, This is both physically and mathematically conceivable Consider the non- 
linear problem, 


V A'VT = pCp-Q^ (3 20) 

where K is taken as a function of temperature The physically conceivable Neumann 
boundary condition to be matched at the interface would be not the derivatives of T but 

0/71 

the fluxes at the interface, — K(T )^- In the case of a linear problem, this boils down to 
the temperature derivatives alone, since K is not a function of temperature Thus for the 
non-lmear problem, global convergence requires 



(3 21) 
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It should be noted that during iterations T\ ^ T<i Also in inhomogeneous media 
with sharp change m thermal conductivity at the interface, heat fluxes rather than the 
temperature gradients should be matched 

3.6.4 Simulation of EOR Using Domain Decomposition 

Let us now consider the simulation of the oil recovery problem for the thm porous formation 
flow domain (Figure 2 1 (a)) using domain decomposition The initial and boundary 
conditions are identical to those considered m Section 2 4 Uzawa’s algorithm is used 
for interface treatment The computer code is thus suitable for parallelization The same 
implicit pressure based formulation as considered for the full domain simulation, is used for 
each sub-domam essentially to circumvent the time step limitations The physical domain 
is split into two equal sub-dimams as shown in Figure 3 4 It is worth mentioning at this 
stage that Uzawa’s algorithm can accomodate any number of sub-domains without any 
additional complications However, we restrict ourselves to two-subdomam formulation to 
make the description mathematically more conceivable 

y 



Figure 3 4 Thm Porous Formation Flow Domain and its Division into Two Subdomains 
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The interface, C, is a stiaighl lme parallel to the y-axis This choice of interface saves 

a considerable amount of effort that will be required while calculating the denvatives of 

the dependent variables along the outward drawn normal Thus for the sub-domains 1 

and 2, boundary conditions on three edges are automatically prescribed from the original 

problem and only the interface has to be taken care of We solve the pressure equations 

for the phase pressures, p 0 and p w , m each sub-domain, also to be specially mentioned 

that the problem is highly non-hneai Hence the quantities to be matched at the interface 

are the phase pressuies, p a and p w and the Darcy velocities of oil and water This is 

in analogy with flux matching as described in Section 3 6 3 The Darcy velocities are, 
'dpo\ — j ( kk r ,A f d% 


kk, 


ZSL 


and 


~^rj Since K and p are constants with respect to the 


Ho ) \ dx ) “““ \p w 
phase pressures of oil and water, considered as Newtonian fluids, we are essentially left 

%■) (tf 


with pressure fluxes, — k r 


, to be matched at the interface De- 


rivatives with respect to x are calculated along the outward normal at the interface for 
each sub-domam A brief outline of the algorithm used m the present study is furnished 
below 

(1) At the beginning of each time step, we supply the guess values of the pressure fluxes 
at the interface The most justified choice would be the fluxes calculated m the previous 
time step Let us call them A* and Aj 

(2) With the chosen values of AJ and A*, we solve the pressure equations m sub-domams 
1 and 2 using AJ and A^ as the boundary conditions at the interface Let the solutions be 
ph and pl t m sub-domam 1 and p l o2 and p x w2 in sub-domam 2 

(3) Using Uzawa’s algorithm, we update the fluxes as 


A o = X l (pLIo -Poilc) 


w (Pw2\c - Pwllc) 


where, to is the interface convergence parameter 
(4) Steps 2-3 are repeated until the conditions 

< S 


|y« - k 
K +i - K 


< s 


are satisfied Here, I indicates the level of interface iteration and 8 indicates the con- 



3 7 Smoothemng of Initial Conditions 


32 


veigence limit for mteiface iterations 

(5) Using the pressure values calculated m sub-domams 1 and 2, S w is updated and ppv 
is calculated 

(6) Fiesh computation is then started for the next time step and the sequence of opera- 
tions, steps 1-5, is repeated 


3.7 Smoothening of Initial Conditions 

The initial conditions for p 0 , p w , S w and T have been smoothened out to avoid compu- 
tational instabilities resulting from a step like initial condition However, these spurious 
oscillations aie observed for short times when initial conditions are not smooth Ini tial 
conditions are of the form given below 

For the thm porous formation domain of nominal size (10 m x 1 2 m) and grid size 
(101 x 11), the initial conditions are as follows 


Pw(x,y) 

= (260 0 

- 70 x) psi 

0 

< X 

< 

0 

1 

X 


= 190 



0 

1 X 

< 

X 

< 

X 

S w (x,y) 

= (0 86 - 

0 66 

x) psi 

0 

< X 

< 

0 

1 

X 


= 0 2 



0 

1 X 

< 

X 

< 

X 

T(x,y) 

= (100 0 

- 125 

x) psi 

0 

< X 

< 

0 

1 

X 


= 50 °C 



0 

1 X 

< 

X 

< 

X 


The initial conditions, for the thick porous formation (the two-spot model) domain 
with a nominal size of (6 m x 6 m) and grid size of (61 x 61), for isothermal injection 
only can be summed up as below 

Pw(x,y) = (260 0 - 116 6 x) psi 0 < x < 0 1 X 
= 190 0 1 X < x < X 


S w {x,y) 


(0 86 - 1 1 x) psi 

0 2 


0 < x < 0 1 X 
0 1 X < x < X 
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In this work, we have implemented the domain decomposition technique only for 
isothermal injection and for the thin porous formation domain of nominal size (4 m x 0 8 
m) and grid size of (81 x 5) However, the initial conditions for the same are as follows 


Pw{x,y) 


(260 0 - 175 x) psi 0 < x < 0 1 X 
190 0 1 X<i<X 


S w (x,y) 


(0 86 - 1 65 x) psi 0<i<01X 
0 2 0 1 X < x < X 


3.8 Code Validation 


The correctness of the computer code has been tested by solving linear and non-linear 
transient heat conduction equations The matrix obtained from the pressure equations has 
been solved by both the specialized sparse matrix solver based on Gaussian elimination and 
by preconditioned conjugate gradient method and the results are seen to be matching quite 
perfectly Also the code has been tested for various grids The domain decomposition code 
has been written m such a manner that it can be parallelized without major modifications 
Calculations in sub-domains are performed m separate subroutines and a mam program 
invokes them and performs the necessary operations The main program has been written 
m a versatile manner so that it can accomodate as many sub-domains as possible 


3.9 Typical Computer Requirements 


In the present study all the runs have been taken on a Dec- Alpha 255/233 workstation 
The isothermal version of the code took around 1 hour and 20 minutes of CPU time for 
a real time simulation of 3 hours for the thin porous formation flow domain and a time 
step (At) of 0 01 hour For the same flow domain and for three hour simulation, the 
non-isothermal version of the code took approximately 1 hour of CPU time for a time step 
(At) of 0 005 hour The isothermal version of the code for the two-spot problem took 
around 2 hours of CPU time for a simulation of three hours and a time step At of 0 01 
hours The convergence criteria used was (Q p+1 — Q p )/Q p+1 x 100 < 0 01, where, Q is 
p 0 , p w and T for all the above cases 



Chapter 4 


Results and Discussions 


In this chapter the results have been grouped under the following sections 

(1) Results related to the flow domain of thin porous formation for both isothermal and 
non-isothermal injection and the effect of surfactants for isothermal injection, 

(2) Various numerical issues corresponding to the domain decomposition technique applied 
to the flow domain of thm porous formation, and 

(3) Results in correspondance with the thick porous formation flow domain i e for the 
two-spot model 


4.1 Results for the Thin Porous Formation Flow 
Domain 

The following issues have been discussed m this section 

1 Effect of formation temperature (T/) on oil displacement efficiency (ppv) 

2 Comparison of oil recovery for the isothermal and non-isothermal cases for a given 
value of Tj 

3 Evolution of saturation and temperature profiles for the isothermal and non-isothermal 
cases 

4 Effect of heat loss, to the overhang and underhang, on oil recovery 
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5 Effect of surfactants on oil recovery 

4.1.1 Isothermal Injection 

Figure 4 1 shows the rise m saturation level of water m the flow domain along its centreline 
with time There is no front and the profile is diffusion dominated Next, we consider the 
effect of increasing the formation temperature We find an increase m the water saruration 
along the centreline of the flow domain as we raise the formation temperature from 30 C 
to 50 °C m steps of ten (Figure 4 2) The improvement in oil recovery with the increase m 
formation temperature has been depicted m Figure 4 3 This is however expected because 
with the increase m temperature, the oil viscosity decreases, thus resulting in an icrease 
m the mobility of oil Figure 4 4 shows the pressure profiles for both oil and water 



Figure 4 1 Progress of Saturation Front for Isothermal Injection of Water 
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Figure 4 2 Effect of Formation Temperature (T/) on Saturation Front for Isothermal 
Injection of Water 



Figure 4 3 Effect of Formation Temperature (T/) on Oil Displacement Efficiency for 
Isothermal Injection of Water 
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Figuie 4 4 Progress of Pressure Fronts of Oil and Water for Isothermal Injection of Water 

4.1 2 Oil Displacement Efficiency 

Oil displacement efficiencies for the isothermal (cold water injection) and non-isothermal 
(hot watei injection) cases are compared m Figure 4 5 The ppv m the non-isothermal 
case is seen to be below that of the isothermal case A possible explanation to this is 
that the accumulation of oil ahead of the water saturation front m the non-isothermal case 
mcieases the oveiall formation resistance However, this phenomenon is not present m 
the isotheimal case The phenomenon leads to an overall decrease m the oil displacement 
efficiency for non-isothermal injection as compared to the isothermal injection, although 
the revei se may be expected elsewhere 
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Figure 4 5 Oil Displacement Efficiency as a Function of Time for Isothermal and Non- 
lsothermal Injection of Water 


4.1.3 Non-isothermal Injection 

The results due to the hot water injection technique are discussed m this subsection 
Figures 4 6 and 4 7 show the development of centreline saturation and temperature profiles 
respectively Both of them reveal a distinct front like structure and on comparing them, we 
see that the saturation front moves more or less with the same speed as the temperature 
front The saturation profile shows an unusual dip immediately beyond the front A 
possible explanation may be that in the vicinity close to the injection zone, due to the 
temperature front the viscous resistance decreases appreciably and the capillary resistance 
is also not very high It essentially increases the waterflood. But m the downstream region, 
both viscous resistance and capillary resistance predominate and an appreciable drop m 
the water saturation values is found Figure 4 8 shows the progress of pressure fronts for 
oil and water phase pressures The oil and water pressures drop steadily after a distance 
of 1 2 m from the injection well 
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Figure 4 8 Progress of Pressure Fronts of Oil and Water for Hot Water Injection 

4.1.4 Effect of Heat Loss 

The effect of heat loss to the overhang and the underhang on oil recovery has been discussed 
m this subsection The pertinent parameter which helps m determining the heat loss is 
Biot number, Bi The higher the value of Bi, the greater is the heat loss Since convection 
plays a dominant role in the transportation of oil, the effect of heat loss does not influence 
the centreline characteristics significantly for the given domain width (Y = 1 2 m) and 
for the time-frame considered here Figure 4 9 shows the typical widthwise temperature 
profiles for various values of Bi For Bi = 0 0 and 0 01, the profiles are straight lines, 
which implies that the overhang and the underhang are perfect insulators The profile for 
Bi = 1 is more or less a straightline and shows up quite significant insulation effect But 
for increasing values of Bi from 100 to 10000, the profile becomes increasingly parabolic 
in nature The slopes show maxima at the top and the bottom of the domain, because of 
the heat loss to the host rocks This effect is dominant near the region very close to the 
injection zone However, m further downstream the effect has been seen to be significantly 
small or not found at all In the figure, the results have been given for a distance of X = 
1 25 m from the injection zone 
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Figure 4 9 Effect of Bi on Widthwise Temperature Profile for Hot Water Injection 

4.1.5 Effect of Surfactants 

The effect of surfactant flooding on oil recovery has been explained in this subsection 
Figure 4 10 shows the oil recovery in terms of ppv as a function of the formation tem- 
perature The two cases, namely, injection without surfactants and with surfactants have 
been compared here It is seen that ppv increases with temperature This is because an 
increase m formation temperature lowers the oil viscosity and hence enhances the mobility 
of oil-water bank The decrease m surface tension in the presence of surfactants is also 
clear in the figure As a consequence, the increase in oil recovery m the presence of surfact- 
ants is quite appreciable as compared to that without surfactants Figure 4 11 compares 
water saturation profiles after three hours of injection with and without surfactants It is 
seen that both the profiles are influenced by the diffusion However, the saturation curve 
for injection with surfactants always shows higher values than that for injection without 
surfactants Thus, a greater mobility of the oil-water bank m the presence of surfactants 
is ascertained Figure 4 12 shows the pressure profiles for oil and water for isothermal 
injection with surfactants The pressure drop for both the components is gradual 
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Figure 4 10 Oil Displacement Efficiency as a Function of Formation Temperature (Tf) 
for Isotheimal Injection of Water with and without Surfactants 



Figure 4 11 Variation of Water Saturation with Distance at the End of Three Hours for 
Isotheimal Injection of Water with and without Surfactants 
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Figure 4 12 Progress of Pressure Fronts of Oil and Water for Isothermal Injection of 
Water with Surfactants 


4.1.6 Application of PCG Method 


The preconditioned conjugate gradient method (PCG) has been applied to solve the matrix 
arising out of the pressure equations for the hot water injection case Figures 4 13 and 4 14 
show the comparison of the saturation and temperature fronts, predicted by the PCG and 
by the Gaussian Elimination, respectively The results are seen to be matching almost 
perfectly This comparison serves to establish the accuracy of the individual solution 
techniques 
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4.2 Numerical Issues Regarding Domain 
Decomposition 


Various numerical issues related to the domain decomposition technique applicable to the 
simulation of flow through the thm porous formation have been discussed herein 

Table 4 1 compares the saturation values calculated using the full domain simulation 
and the simulation with domain decomposition The values are seen to be sufficiently close 
at various distances from the injection zone Table 4 2 shows the comparison between 
the CPU times taken foi five hour simulation using the full domain simulation and the 
domain decomposition code Domain decomposition is seen to take more CPU time 
This can be attributed to the use of a sequential machine m solving a highly non-lmear 
problem Implementation of parallelization is expected to show an overall speed up of the 
computing time Also to be mentioned that for the linear problems, domain decomposition 
has been found to be faster even on sequential machine For both the cases, the formation 
temperature is 30° C and the domian decomposition algorithm has been applied for two 
sub-domains 


Table 4 1. Comparison of Saturation Values for Full Domain Simulation and Simulation 
with Domain Decomposition 



x, m 

05 

1 0 

1 5 

20 

25 

30 

35 

40 


A* 

0 3595 

0 2837 

0 2589 

0 2538 

0 2016 

0 2009 

0 2113 

0 2000 

II 

B f 

0.3513 

0 2789 

0 2617 

0 2420 

0.2717 

0 2451 

0 2175 

0 2004 

t — 5 hr 

A* 

0.4585 

0 3379 

0 3287 

0 2872 

0 2637 

0 2404 

0 2585 

0 2000 


B T 

0 4500 

0 3355 

0 3214 

0 2916 

0 2578 

0 2300 
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0 2004 


* Full domain simulation 
f Simulation with domain decomposition 
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Table 4 2 Comparison of CPU 1 mit \ alues for Full Domain Simulation and Simulation 
with Domain Decomposition 



CPU time 

A 

130 seconds 

B 

3560 seconds 


* Full domain simulation 
f Simulation with domain decomposition 


4.2.1 Effect of Interface Convergence Parameter 

Here we are discussing the effect of interface convergence parameter (a;) on domain de- 
composition We know that u is an adjustable constant that lies between 0 and 1 In 
the present work we have performed a thorough numerical experimentation to find out 
an optimum value of u However, the subsequent text shows a comparative study for 
three values of the interface parameter (u> = 0 2, 0 25 and 0 3) based on the CPU time 
taken for a five hour simulation (Table 4.3) Also, the ppv values at the end of five hours 
with respect to different u> values have been compared with the results due to full domain 
simulation (Table 4 4) 


Table 4 3. Comparison of CPU Time Values 



CPU time 

u = 0 2 

3800 seconds 

u = 0 25 

3560 seconds 

u> = 0 3 

3900 seconds 

"FDS* 

130 seconds 


| Full Domain Simulation 
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Table 4 4 Comparison of ppv Values 


t, 

hr 



i 


2 


3 

4 


5 

u 

= 0 

2 

2 

5447 

4 

9307 

7 

8495 

10.016 

12 

739 

U) 

= 0 

25 

2 

8382 

5 

0785 

7 

2961 

10 168 

11 

911 

UJ 

= 0 

3 

2 

8369 

5 

1859 

7 

8202 

10.110 

12, 

.083 

FDS* 

2 

8299 

5 

3882 

7 

7988 

10 175 

12, 

.536 


| Full Domain Simulation 


Fiom the above comparison, it is clear that the ppv values for cu = 0 25 and 0.3 axe 
very close to that obtained fiom the full domain simulation. But the CPU time for u = 
0 25 is less than that for lo = 0 3. For the problem of interest, we have identified that w = 
0 25 is the optimum value and all the calculations concerning the domain decomposition 
technique have been done based on that value. 


4.2.2 Effect of Uzawa Convergence Limit 

The effect of convergence limit for the interface iterations on the simulation with domain 
decomposition as compared to the full domain simulation has been discussed in this sub- 
section Figures 4 15 and 4.16 show the pressure profiles for oil and water respectively for 
two values (10~ 3 and 10~ 4 ) of the Uzawa Convergence Limit (hereafter UCL) and those 
for the full domain simulation For both the cases, the interface convergence parameter 
(w) is 0 25 and the formation temperature is 30° C. The figures clearly show that for a 
lower value of the UCL, the pressure profiles are closer to those obtained from the full 
domain simulation. It will obviously predict a much better ppv value at the end of five 
hour simulation as it is observed from the Table 4.5. Table 4 5 shows a comparison of the 
ppv values for simulation with domain decomposition for the two UCL values and the full 
domain simulation 
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Figure 4 15 Comparison of Oil Pressure Profile for Various Values of Uzawa Convergence 
Limit and for Full Domain Simulation 



x, m 


Figure 4 16 Comparison of Oil Pressure Profile for Various Values of Uzawa Convergence 
Limit and for Full Domain Simulation 
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Table 4 5. Comparison of ppv Values 



ppv 

UCL* = IQ” 3 

11.911 

UCL = IQ- 4 

T2 208 

FDS* 

12 536 


* Uzavja Convergence Limit 
• Full Domain Simulation 


4.2.3 Effect of Interface Treatment 

This subsection discusses the effect of various interface treatments on simulation with 
domain decomposition. In the present study, we have considered a variant of three inter- 
face treatment approaches In the first approach, we have taken only one-sided normal 
derivative of the phase pressure gradients at the interface In the second approach, we 
have considered the averaging of the normal derivatives of the phase pressure gradients 
at the interface However, both the approaches give very close ppv values for a five hour 
simulation, as shown in Table 4.6. In both the cases, the phase pressure gradients are 
used as the initial guess for the interface condition and are updated subsequently through 
Uzawa iterations. Therefore, the final solution does not change In the third approach, we 
have used the mass fluxes, derived from the Darcy’s law (as discussed in Section 3.6.4), in 
stead of the pressure gradients However, due to some unidentified difficulties experienced 
m implementing the mass flux approach, we are unable to furnish any result here. The 
failure forms a strong motivation for the future work using the mass-flux based interface- 
treatment. In the first and second approaches, we have taken u> = 0.25, Tj = 30°C and 
UCL = 10" 3 


Table 4 6: Comparison of ppv Values 



ppv 

FA* 

11 911 

“SA^ 

12.022 


o First Approach 
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4.3 Results for the Two-spot Model 


In this section, various issues t elated to the simulation of the two spot model have been 
discussed Results foi the thieu simplified versions (Figures 2 2 (ck f d > and (e)) of the 
original two-spot model (2 2 (hi) have been piesented here 


4.3.1 Oil Displacement Efficiency 

Figures 4 17, 4 18 and 1 19 show the oil displacement efficiencies for the isothermal in- 
jection with and without surfactants The ppv values are seen to be higher for injection 
with surfactants This is expected as surfactants reduce the surface tension effect, thus 
enhancing the mobility of the oil-water bank Two values of the formation temperature, T f 
= 40° C and 50°C, have been considered Figures 4 17 (a), 4 18 (a) and 4 19 (a) show the 
displacement efficiencies for 7/ = 40 o (’ and Figures 4 17 (b), 4.18 (b) and 4.19 (b) show 
the displacement efficiencies for 7/ = 5Q°C. However, the figures show that for surfactant 
flooding (for the cases corresponding to Figure 2.2 (c) and Figure 2 2 (e)), it was not 
possible to go for full five hour simulation. A possible explanation behind this peculiar 
phenomenon is that surfactants enhance the mobility of the oil-water bank to such an 
extent that from some regions in the flow domain oil is fully flushed out and only water 
remains This implies that after a few hours of simulation, in some regions the water 
saturation, S w becomes 1. In other words, oil saturation, S 0 tends to 0. Our model is 
not valid for such a situation. We assume that the irreducible limit for oil is S w = 0 86 
and that for water is S w = 0.1 Therefore, for the simulation of surfactant flooding, we 
require to choose the formation temperature in such a way that the irreducible limit is not 
overlooked It may be mentioned that full five hour simulation was accomplished for the 
case corresponding to Figure 2.2 (d) and the formation temperature, Tj = 40° C (see the 
results in Figure 4.18 (a)) 

4.3.2 Transverse Pressure Profile 

Figures 4 20, 4.21 and 4,22 show the evolution of the transverse pressure profiles for the 
flow domains as shown in Figures 2.2 (c), (d) and (e) respectively. The profiles are seen 
to exhibit a gradual increasing trend in slope along the y direction. This is, however, 
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expected for the sustenance of the flow of oil from the injection well to the production 
well For all the above three cases the formation temperature has been taken as 50°C and 
isothermal injection without surfactants has been considered Figures 4 20 (a), 4 21 (a) 
and 4 22 (a) show the pressure profiles for the oil phase and Figures 4 20 (b), 4 21 (b) and 
4 22 (b) show the pressure profiles for the water phase 

4.3.3 Saturation Profile 

The distribution of saturation (S w ) m the flow domain for the two-spot model has been 
predicted here Figure 4 23 shows the saturation distribution for the flow domain (cor- 
responding to Figure 2 2 (c)) for various durations of simulation Figure 4 23 (a) shows 
the initial distribution of S w while Figures 4 23 (b), 4 23 (c) and 4 23 (d) predict the 
distribution of S w after three hours, five hours and ten hours respectively For the other 
two flow domains (as shown in Figures 2 2 (d) and (e)) more or less same trend of vari- 
ation of the saturation distribution has been observed However, from the figures, it can 
be conjectured that S w increases unusually in the region of the production well for the 
simulations for three hours, five hours and ten hours This may be due to the intrinsic 
approximation of the original two-spot model which entails accumulation of water near 
the impermeable bounding surfaces. However, the accumulation of water m the vicinity 
of the production well in the present simulation has been unusually high Also here, we 
have taken the injection and production well diameters as 1 meter while the dimension 
of the two-spot model flow domain (Figure 2 2 (b)) is only 6 meters But in real life, 
the two-spot domain size itself may be several times of the domain size considered in the 
present study This may also lead to unusual behaviour in the saturation distribution 
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(a) 



Figure 4 17 Oil Displacement Efficiency as a Function of Time for Isothermal Injectioi 
of Water with and without Surfactants for the Two-spot Model Flow Domain as shown n 
Figure 2 2 (c) 
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Figure 4 18 Oil Displacement Efficiency as a Function of Time for Isothermal Injection 

of Water with and without Surfactants for the Two-spot Model Flow Domain as shown in 
Figure 2 2 (d) 
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With Surfactants 
Without Surfactants 


With Surfactants 
Without Surfactants 


Figure 4 19 Oil Displacement Efficiency as a Function of Time for Isothermal Injection 
of Water with and without Surfactants for the Two-spot Model Flow Domain as shown m 
Figure 2 2 (e) 
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Figure 4 20 Variation of Transverse Pressure Profiles of Oil and Water for Isothermal 
Injection of Water without Surfactants for the Two-spot Model Flow Domain as shown m 
Figure 2 2 (c) 
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(a) 



Pressure, MPa 
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Figure 4 21. Variation of Transverse Pressure Profiles of Oil and Water for Isothermal 
Injection of Water without Surfactants for the Two-spot Model Flow Domain as shown m 
Figure 2 2 (d) 
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Figure 4 22 Variation of Transverse Pressure Profiles of Oil and Water for Isothermal 
Injection of Water without Surfactants for the Two-spot Model Flow Domain as shown in 
Figure 2 2 (e) 
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4.4 Concluding Remarks 


Before concluding, we would like to mention about role of the field data in the form 
of constitutive relations as furnished in Table A 1 of Appendix A In the table, k ro w is 
given for a finite set of S w values A linear interpolation is used to compute intermediate 
values of the relative permeability and capillary pressure Also it is well known that the 
prediction of a numerical simulator for flow through a porous media is sensitive to the 
values of k rcn k rw and p Cow . Any deviation from the correct values with respect to these 
input parameters can lead to the unrealistic predictions. However, the pointers from the 
open literature suggest that the closed form analytical expressions are expected to improve 
the quality of the numerical predictions (Pruess and Narasimhan, 1982) 

Also it is essential to point out the relative advantages of harmonic averaging over 
arithmetic averaging of 6 and 7 values as found in the pressure equations (Section 3.1) 
The reason behind the choice of harmonic averaging is that the viscosity has very sensitive 
exponential dependence on temperature and the relative permeability changes sharply with 
saturation The resultant effect can cause sharp variations in 6 and 7 values in a small 
region However, the choice proved to be a correct one as the wiggles m saturation profiles 
for the isothermal case were remarkably reduced in magnitude when the property values 
were determined using harmonic averaging procedure. 


Chapter 5 


Conclusions and Scope for Future Work 


In this work, we have tried to take up the numerical simulation of the enhanced oil recovery 
problem from different perspectives. It includes modeling and simulation of flow through 
two different domains The two-spot problem is a major venture in this study It gives 
an idea as to how field-scale problems can be tackled The implementation of the domain 
decomposition technique is yet another attempt to study the oil recovery problem from 
the stand point of a real reservoir simulation. It forms the building block of simulation of 
the large reservoirs The implementation of the PCG method m solving the matrix aris in g 
out of the pressure equations is another salient feature of the present work and it is a 
primafacie attempt towards parallelizing the complete code in future 

Recommendations 

The following extensions to the present work are recommended to make the numer- 
ical simulation of the oil recovery problem complete, robust and realistic 

(1) Stability limits of simulation: 

The numerical work presented in this work has been done for one particular forma- 
tion i e one particular set of formation and surfactant properties Due to high degree of 
non-linearity in the governing equations, application of this simulation code to other sets of 
properties cannot be guaranteed. Hence to ascertain versatahty of the present formulation 
and computer code, simulation should be performed with other variants of the formation 
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(2) Three-phase flow: 

In addition to oil and water, the third phase present in certain reservoirs is air, 
though comparatively much less m quantity A three phase simulation of oil, water and 
air would thus lend more credibility to this work Also we can take up the three phase 
simulation of oil, water and steam for the future analysis 

(3) Adsorption-Desorption in the formation: 

i 

In the case of surfactant flooding, we can take into account the adsorption-desorption 
phenomena, i.e transfer of surfactants from the injected flow to the formation and vice 
versa respevtively The mass transfer phenomena will obviously call for modification in the 
constitutive relations in conjunction with the variation of concentration of the surfactant 
solution It will lead to the inclusion of an additional mass transfer equation that must be 
solved along with the pressure equations 

(4) Enhancement in the two-spot problem: 

In case of the two-spot problem, we can consider the original two-spot model (Figure 
2 2 (b)), m stead of the simplified flow domains (Figures 2 2 (c), (d) and (e)) Obviously, 
it will include variable grid and proper handling of the boundary conditions on the circular 
axe shaped injection and production zones. It will invariably make the problem compu- 
tationally more intensive Enhancement of this problem for non-isothermal case is also 
highly recommended for further analysis. It will invoke an additional energy equation 
which must be solved along with the pressure equations, thus making the problem much 
more computation expensive Also we know that in the vicinity of the injection and pro- 
duction wells, near-field effects are predominant. Therefore, we can extend the domain 
decomposition technique by localizing the affected regions through assigning separate sub- 
domains to them and studying them closely 

(5) Three dimensional simulation: 

An extension of the present work for the rectangular flow domain can be a three 
dimensional simulation including mhomogeneity and anisotropy of the formation This 
will also make the reservoir simulation computationally more expensive. Domain decom- 
position with parallelization can be a very good option for this problem. 
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(6) Interface treatment: 

In the domain decomposition problem, Uzawa’s algorithm should be looked into 
much more details for the interface treatment In this regard, matching of the mass fluxes 
instead of the pressure gradients at the interface, is highly recommended Also we can 
go for some advanced parallelizable techniques ( Glowinski, 1983), in lieu of the Uzawa’s 
algorithm which can be expectecd to enhance the convergence rate of interface iterations 
and make the entire simulation faster 

(7) Parallelization: 

The major task recommended for future extension is to parallelize the code developed 
for the domain decomposition problem. The successful implementation would not only save 
a lot of computational time but also enable us to perform real life reservoir simulation 
problems. 
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Table A.l Fluid and Formation Properties 


K 

132 darcies 

Kh 

0 1661 W/m °C 

e 

0 375 

£<> 

0 03447 Pa" 1 


0 02137 Pa" 1 

Po 

2 28 x lO" 4 °C- 1 

Pu, 

2 28 x 10~ 4 °C -1 

C 0 

2092 (J/kg-°C) 

c w 

4184 (J/kg-°C) 

(p c )r 

2413 (J/m 3 -°C) 

Po 

90% of p w 


Table A 2 Oil and Water Viscosities 


o 

O 

20 

40 

60 

80 

Po, Pa-s 
p w , Pa-s 

4.12 

1.03 x 10" 4 

2.37 

6.6 x 10" 4 

0.61 

4 73 x 10 -4 

0.10 

3.67 x lO" 4 


Table A.3. Constitutive Relations 


Sw 

k 

f^rw 

Ko 

p Cew , KPa 

0 1* 

0 

1.0 

28.34 

02 • 

0.0016 

0 875 

0 655 

0.3 

0 0081 

0.735 

0 496 

04 

0 0259 

0 590 

0 4205 

0.5 

0 0672 

0.420 

0.352 

06 

0.1000 

0.210 

0 283 

0.7 

0 1400 

0 070 

0 214 

08 

0 2000 

0.016 

0.145 

0 86+ 

0 2500 

0 

0 0758 


* Irreducible limit for water 
f Irreducible limit for oil 



